#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#include<time.h>
#define PI 3.1415926
#define N 1024
#define SHOW_N 32

typedef struct {
    float re;// really
    float im;// imaginary
} complex,*pcomplex;

static void func_write_txt(const char *file_name, const char *write_mode, int sample_count, float *data);

complex complexBuild(float re,float im)
{
    complex cx;
    cx.re=re;
    cx.im=im;
    return cx;
}


complex complexAdd(complex a,complex b)
{
    complex ret;
    ret.re=a.re+b.re;
    ret.im=a.im+b.im;
    return ret;
}

complex complexMult(complex a,complex b) {
    complex ret;
    ret.re=a.re*b.re-a.im*b.im;
    ret.im=a.im*b.re+a.re*b.im;
    return ret;
}
void DFT(complex x[],complex X[], float *dft_amp)
{
    int k,n;
    complex Wnk;
    for (k=0; k<N; k++)
    {
        X[k].re=0;
        X[k].im=0;
        for (n=0; n<N; n++)
        {
            #if 0
            Wnk.re=(float)cos(2*PI*k*n/N);
            Wnk.im=(float)-sin(2*PI*k*n/N);
            X[k]=complexAdd(X[k],complexMult(x[n],Wnk));
            #else 
            //Wnk.re=(float)cos((-2)*PI*k*n/N);
            //Wnk.im=(float)sin((-2)*PI*k*n/N);
            //X[k]=complexAdd(X[k],complexMult(x[n],Wnk));
            //complex tmp;
            //tmp = complexMult(x[n],Wnk);
            //tmp.re = x[n].re * Wnk.re - x[n].im * Wnk.im; //complexMult(,Wnk);
            //tmp.im = x[n].im * Wnk.re + x[n].re * Wnk.im; // complexMult(x[n],Wnk);
            X[k].re = X[k].re + x[n].re * cos((-2)*PI*k*n/N) - x[n].im * sin((-2)*PI*k*n/N);
            X[k].im = X[k].im + x[n].im * cos((-2)*PI*k*n/N) + x[n].re * sin((-2)*PI*k*n/N);
            #endif
        }
        dft_amp[k] = sqrt(pow(X[k].re, 2) + pow(X[k].im, 2));
    }

}

void IDFT(complex X[],complex x[])
{   int k,n;
    float im=0;
    complex ejw;
    for (k=0; k<N; k++)
    {
        x[k].re=0;
        x[k].im=0;
        for (n=0; n<N; n++)
        {
            #if 0
            ejw.re=(float)cos(2*PI*k*n/N);
            ejw.im=(float)sin(2*PI*k*n/N);
            x[k]=complexAdd(x[k],complexMult(X[n],ejw));
            #else
            x[k].re = x[k].re + X[n].re * cos((-2)*PI*k*n/N) + X[n].im * sin((-2)*PI*k*n/N);
            x[k].im = x[k].im + X[n].im * cos((-2)*PI*k*n/N) + X[n].re * sin((-2)*PI*k*n/N);
            #endif
        }
        x[k].re/=N;
        x[k].im/=N;
    }
}
static void func_write_txt(const char *file_name, const char *write_mode, int sample_count, float *dataR)
{
    FILE *fp = NULL;
    int i= 0;
    if((fp=fopen(file_name, write_mode))==NULL)
    {
        printf("[%s +%d %s] fopen %s failed!\n", __FILE__, __LINE__, __func__, file_name);
        exit(0);
    }
    for(i = 0; i < sample_count; i++)
        fprintf(fp,"%f,",dataR[i]);
    fclose(fp);
}
void show_complex_messages3(float *real, float *imag, float *amp)
{
    for(int i = 0; i < SHOW_N; i++)
    {
        printf(" (%f, %f ) %f ", real[i], imag[i], amp[i]);
    }
    printf("\n");
}
void show_complex_messages2(float *real, float *imag)
{
    for(int i = 0; i < SHOW_N; i++)
    {
        printf(" (%f, %f) ", real[i], imag[i]);
    }
    printf("\n");
}
void show_complex_messages1(complex *input)
{
    for(int i = 0; i < SHOW_N; i++)
    {
        printf(" (%f,%f) ",input[i].re,input[i].im);
    }
    printf("\n");
}
void show_float_messages3(float *real, float *imag, float *amp)
{
    for(int i = 0; i < SHOW_N; i++)
    {
        printf(" (%f, %f ) %f ", real[i], imag[i], amp[i]);
    }
    printf("\n");
}
void show_float_messages2(float *real, float *imag)
{
    for(int i = 0; i < SHOW_N; i++)
    {
        printf(" (%f, %f) ", real[i], imag[i]);
    }
    printf("\n");
}
void show_float_messages1(float *input)
{
    for(int i = 0; i < SHOW_N; i++)
    {
        printf(" (%f) ", input[i]);
    }
    printf("\n");
}
int main()
{
    complex samples[N],_out[N];
    float input[N];
    float dft_amp[N];
    float idft[N];
    int i;
    int k=1;
    int f1 = 100;//   	# 信号的特征频率1
    int f2 = 300;//   # 信号的特征频率2
    int f3 = 500;//   # 信号的特征频率3
    for (i=0; i<N; i++)
    {
        samples[i].re=i;
        samples[i].im=0;
        float j = i;
        samples[i].re = 10 + 7 * sin(2*PI*f1*(j*1/N)) + 5 * sin(2*PI*f2* (j*1/N)) + 3 * sin(2*PI*f3* (j*1/N));
    }
    printf("cpmplex samples:\n");
    show_complex_messages1(samples);
    for(int i = 0; i < N; i++)
    {
       input[i] = samples[i].re;
    }
    printf("input samples:\n");
    show_float_messages1(input);
    func_write_txt("complex_dft_input.txt", "wb", N, input);

    DFT(samples,_out, dft_amp);
    printf("comlex DFT(%d):\n", N);
    show_complex_messages1(_out);

    printf("complex_dft_amp(%d):\n", N);
    show_float_messages1(dft_amp);

    func_write_txt("complex_dft_output.txt", "wb", N, dft_amp);
    for(int i = 0; i < SHOW_N; i++)
    {
        printf(" (%f) ",dft_amp[i]);
    }
    printf("\n");
    IDFT(_out,samples);
    printf("comlex IDFT(%d):\n", N);
    show_complex_messages1(samples);
    for(int i = 0; i < N; i++)
    {
       idft[i] = samples[i].re;
    }
    func_write_txt("complex_idft_output.txt", "wb", N, idft);
    return 0;
}
